##=============================================================================
## Appendix Figure 18
##=============================================================================

##-----------------
# clear environment
rm(list=ls())
options(stringsAsFactors = FALSE, scipen = 999)
# source("R/functions.R")

seed <- sample.int(.Machine$integer.max, 1)
set.seed(seed)

ipak <- function(pkg){new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "ggpubr")

ipak(packages)

##---------
# Load data
#setwd("/Users/austinknuppe/Library/CloudStorage/Dropbox/Ukraine2022WartimeSurvey/Paper_peace/final_version_oct_2024/replication-scripts")
load("~/Desktop/clean_ukraine_data.RData")

##----------------------------------------
# Sensitivity Analysis with Imputed Values

# https://www.gerkovink.com/miceVignettes/Sensitivity_analysis/Sensitivity_analysis.html

dat <- dat %>% 
  select(id, wave, adm1_en, survival_mindset, save_lives, compromise,
         compromise_idx_full, sex, age, displaced, ethnic, home_lang, edu, job, finances)
  
# Perform a dry run (using maxit = 0) in mice. 
ini <- mice(dat, maxit = 0)

#  Generating imputations under the δ-adjustment
delta <- c(0, -0.05, -0.10, -0.15, -0.20, 0.05, 0.10, 0.15, 0.20)

imp.all <- vector("list", length(delta))

post <- ini$post

# NB replace for-loop with function for better comp. speed
for (i in 1:length(delta)){
  d <- delta[i]
  cmd <- paste("imp[[j]][,i] <- imp[[j]][,i] +", d)
  post["survival_mindset"] <- cmd
  imp <- mice(dat, post = post, maxit = 5, seed = i, print = FALSE)
  imp.all[[i]] <- imp
}

#bwplot(imp.all[[1]]) 
imp_1_plot <- densityplot(imp.all[[1]], lwd = 3) # delta = 0
imp_1_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = 0)"))

#bwplot(imp.all[[5]]) 
imp_5_plot <- densityplot(imp.all[[5]], lwd = 3) # delta = -0.20
imp_5_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = -0.20)"))

#bwplot(imp.all[[9]]) 
imp_9_plot <- densityplot(imp.all[[9]], lwd = 3) # delta = 0.20
imp_9_plot[["xlab"]] <- expression(paste("Survival Mindset (", delta, " = 0.20)"))

ggpubr::ggarrange(imp_1_plot, imp_5_plot, imp_9_plot,
                  ncol = 1, 
                  nrow = 3,
                  labels = c("Imputation 1: Missing Completely at Random", 
                             "Imputation 5: Skewed -1 SD", 
                             "Imputation 9: Skewed +1 SD"))

rm(list = ls())
##=============================================================================
## End of File
##=============================================================================